Phase Field Modeling of Fast Crack Propagation 
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We present a continuum theory which predicts the steady state propagation of cracks. The 
theory overcomes the usual problem of a finite time cusp singularity of the Grinfeld instability by 
the inclusion of elastodynamic effects which restore selection of the steady state tip radius and 
velocity. We developed a phase field model for elastically induced phase transitions; in the limit of 
small or vanishing elastic coefficients in the new phase, fracture can be studied. The simulations 
confirm analytical predictions for fast crack propagation. 

PACS numbers: 62.20.Mk, 46.50.+a, 46.15.-x, 47.54.+r 
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Understanding the day-to-day phenomenon of fracture 
is a major challenge for solid state physics and materi- 
als science. Starting with the early idea of Griffith 0, 
who realized that crack growth is a competition between 
a release of elastic energy and an increase of surface en- 
ergy, various approaches have been developed to describe 
the striking features of cracks 0. Usually, the motion 
of cracks is understood on the level of breaking bonds 
at sharp tips, and obviously theoretical predictions de- 
pend sensitively on the underlying empirical models of 
the atomic properties (see for example |3|). Plastic ef- 
fects, however, lead to extended crack tips (finite tip ra- 
dius r*o), and it is conceivable that for example fracture 
in gels can be described macroscopically. Then a full 
modeling of fracture should not only determine the crack 
speed but also the crack shape self-consistently. 

Recent phase field models go beyond the microscopic 
limit of discrete models with broken rotational symmetry, 
and encompass much of the expected behavior of cracks 
0, ; these models are close in spirit though different in 
details with respect to earlier approaches [fj. However, 
the scale of the appearing patterns is always dictated 
by the phase field interface width, and thus these models 
have problems in the sharp interface limit. Other descrip- 
tions are based on macroscopic equations of motion but 
suffer from inherent finite time singularities which do not 
allow steady state crack growth unless the tip radius is 
limited by the phase field interface width Q. Numerical 
approaches which are not based on a phase field provide a 
selection mechanism by the introduction of complicated 
nonlinear terms in the elastic energy for high strains in 
the tip region ||, requiring additional parameters. 

It is therefore highly desirable to look for minimal mod- 
els of fracture which are free from microscopic details and 
which are based on well established thermodynamical 
concepts. This is also motivated by experimental results 
showing that many features of crack growth are rather 
generic 0; among them is the saturation of the steady 
state velocity appreciably below the Rayleigh speed and 



a tip splitting for high applied tension. 

Already in our previous publication we empha- 
sized a connection between fracture mechanics and elas- 
tically induced surface diffusion processes: the Asaro- 
Tiller-Grinfeld (ATG) instability [HI appears to be a 
good starting point for the quest for a "macroscopic" 
theory of fracture, where crack growth is mediated by 
surface diffusion along the crack surfaces. Plasticity ef- 
fects are contained in this description in a lubrication 
approximation, provided that the region of inelasticity is 
relatively thin in comparison to the tip radius. 

Here we propose a similar approach which describes 
crack growth by a phase transition model. In fact, it 
produces a non-trivial dynamical selection of the radius 
of the crack tip. Instead of cracks filled with vacuum, we 
consider for the moment a soft condensed phase inside 
the crack which is growing at the expense of a brittle 
material pfl. 

The difference in the chemical potentials between two 
phases at an interface is 01 
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provided that the soft phase is stress free because of neg- 
ligible elastic moduli. Then the surface of the crack is 
free of normal and shear stresses. We assume for sim- 
plicity the mass density p to be equal in both phases 
and the elastic displacements to be continuous at the 
interface, which means that the two phases are coher- 
ent. Also, we assume a two-dimensional geometry. The 
interfacial energy per unit area is 7, and the interface 
curvature k is positive if the crack shape is convex. f2 is 
the atomic volume, ajf. and stress and strain tensor 
respectively. Stress and strain are connected by Hooke's 
law for isotropic elasticity, Okj = 2fmkj + Xuudkj, with 
the Lame coefficient A and the shear modulus fi. Alter- 
natively, we use Young's modulus E = p(3\ + 2fi) / (\+ p) 
and the Poisson ratio v = X/2(X+fi) as elastic constants. 
For phase transitions, the motion of the interface is 



locally expressed by the normal velocity 



with a kinetic coefficient D with dimension [D] = m 2 s _1 . 

It is known that nonhydrostatic stresses P at a nom- 
inally flat interface lead to the ATG instability: Long- 
wave perturbations of a flat interface diminish the to- 
tal energy of the system, whereas short-wave perturba- 
tions are hampered by surface energy. The characteristic 
length scale of this instability, Lq ~ E^/P 2 (l — v 2 ), is 
identical to the Griffith length of a crack, up to a numer- 
ical prefactor. This instability leads to a finite time cusp 
singularity in the framework of the static theory of elas- 
ticity [13j: The tip radius decreases to zero and simulta- 
neously the tip velocity grows indefinitely. In this sense, 
the advancing notches can be interpreted as cracks. We 
explained the unphysical breakdown of the late stage of 
the ATG instability in our previous publication for sur- 
face diffusion [lOj: Here, similar arguments apply, and 
detailed counting arguments will be derived below prov- 
ing that simultaneous selection of the tip radius and ve- 
locity is impossible in the framework of the static theory 
of elasticity. 

The elastic problem of a semi-infinite mathematical 
cut in an infinite (two dimensional) solid is exactly 
solved by a square root singularity of stresses, cry = 
K fo,ij{6) I '(27rr) 1 ' 2 , using polar coordinates as depicted 
in Fig. n] Here K is the stress intensity factor (static 
or dynamic), fo,ij the universal angular distribution de- 
pending only on the mode of loading (for brevity, we sup- 
press the velocity dependence v/vr; vr is the Rayleigh 
speed); we concentrate on cracks of "mode I" -type here. 
For an extended crack in an infinite environment, the 
square root singularity is only the slowest decaying mode; 
other powers can be present, and the total field can be 
interpreted as an expansion in the set of eigenfunctions 
of a straight crack: 

(2^)17? E r n • ( 3 ) 
v ' n=0 

Far field conditions imply Aq = 1, whereas the other co- 
efficients A n are determined by the boundary conditions 
of vanishing shear and normal stress on the crack; modes 
with ascending powers of r are forbidden by boundary 
conditions, and even a homogeneous stress P cannot be 
present in an infinite system (divergent strip width L), 
because this would lead to a diverging stress intensity 
factor K ~ PL 1 / 2 . All higher order corrections in Eq. 
P|) scale as A n ~ Tq and vanish in the sharp tip limit. 
From this equation follows readily that the stresses on 
the crack surface scale as a ~ r 1 ^ 2 . 

We consider a crack as it might have developed in the 
late stage of the ATG instability, as depicted in Fig. ^ 
The macroscopic length of the crack is not considered 
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FIG. 1: Steady state growth of a crack in a strip, as obtained 
from phase-field simulations. The grey color corresponds to 
the solid phase. A constant displacement is prescribed at the 
upper and lower boundary of the system. The parameters 
used here are A = 2.03 and D/tVR = 6.18; the system size 
is 600 x 400. The resulting velocity is v/vr — 0.71 and the 
radius ro = 0.69D/vr = 4.26e are dynamically selected. 



here, and instead the stress intensity factor K is given. 
At first, we demonstrate that steady state growth using 
only the static theory of elasticity is impossible; in fact, 
this is the reason for the mentioned cusp singularity. 

Assume that y{x) in Cartesian or r(6) in polar coordi- 
nates describes the steady state shape of a crack in the 
co-moving frame of reference, corresponding to a specific 
tip radius ro and velocity v. According to the results 
above, both contributions to the chemical potential Eq. 
(JTJ scale as fi ~ 1 /ro and thus by virtue of the equation 
of motion J5J), v ~ 1/ro- Hence a rescaling of the steady 
state equation is possible: Increasing the crack size by a 
factor a simply reduces the steady state velocity by 1/a 
and vice versa. In other words, the explicit length scale 
ro drops out of the equations, and only vtq/D remains as 
dimensionless parameter. All other parameters combine 
to the dimensionless driving force A = K 2 (l — v 2 ~)j2E^\ 
A = 1 corresponds to the Griffith point. Notice that this 
rescaling is only possible in the framework of the static 
theory of elasticity; otherwise, the stress field itself be- 
comes velocity dependent, introducing the ratio v/vr as 
additional parameter in the system. 

Using the steady state condition y = —vy', the shape 
equation reads 

OikUik vy' 

K = ~ + D(l + ^"- " 

which is a nonlocal equation due to the long range elas- 
tic interactions. The boundary conditions at the tip are 
given by the arbitrary choice of the origin, r(0) = r , and 
r'(0) = 0, since we are interested in symmetrical shapes, 
r(6) = r(—0). Thus the entire shape is a function de- 
pending only on the parameter v. 

On the other hand, in the tail region, the elastic 
stresses have decayed, and the shape equation there- 
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fore becomes simply —vy' = Dy" . Its general solution, 
y(x) — A + B exp(— ux/D), contains a growing exponen- 
tial which is inconsistent with the boundary conditions 
of the straight crack. Therefore the solution must be ar- 
ranged such that B = 0. Notice that in contrast to the 
surface diffusion process ^3 a finite opening A cannot be 
excluded since we do not have to obey mass conservation 
here. At a first glance, suppression of the exponential 
seems to lead to a selection of the steady state velocity 
v as only available parameter. However, from the rescal- 
ing behavior explored above it follows immediately that 
B ~ 1/v (notice that B has the dimension of a length) , 
and therefore a finite velocity cannot be selected. Con- 
sequently, a steady state solution for a growing crack in 
the framework of the static theory of elasticity does not 
exist. 

The situation is different for the dynamical theory of 
elasticity: Here velocity enters into the equation of mo- 
tion not only as vvq/D but also as v/vr. Hence a second 
independent parameter exists to guarantee B = 0. Since 
now a rescaling is impossible, both the propagation ve- 
locity v and the tip radius ro are selected. 

A tip splitting is possible for high applied tensions due 

— 1/2 

to a secondary ATG instability: Since a ~ Kr Q in 
the tip region and the local ATG length is Lq ~ Ej/a 2 , 
an instability can occur, provided that the tip radius be- 
comes of the order of the ATG length. In dimensionless 
units, this leads to the prediction A sp i it ~ 0(1). 

We developed a phase field code together with elasto- 
dynamics to describe phase transformations under stress, 
including for example also martensitic transformations. 
In the limit of vanishing shear modulus in one of the 
phases, this approach can be used to study melting and 
solidification processes which are induced by elastic forces 
0- For a very soft secondary phase, crack propagation 
can be studied in the framework of a continuum the- 
ory, since then the usual boundary conditions of vanish- 
ing normal and shear stress are recovered. Let <fi de- 
note the phase field with values = for the soft and 
4> = 1 for the hard phase. The energy density contri- 
butions are f e i = ^i(cf))u 2 j + X((f>)(uu) 2 /2 for the elastic 
energy, with fj,(4>) = h(4>)^ + (1 - h(cj)))^ 2 '> and X(cf) = 
h(<f>)\W+(l-h(<f>))\W, where h(<j>) = 2 (3-20) interpo- 
lates between the phases and the superscripts denote the 
bulk values. The surface energy is f s {4>) = 37e(V0) 2 /2 
with the interface width e. Finally, fd w = 6~f(j> 2 (l — 4>) 2 l e 
is the well-known double well potential. Thus the total 
potential energy is 

U = J dV (fel + fs + fdw) ■ (5) 

The elastodynamic equations are derived from the energy 
by the variation with respect to the displacements Ui, 

pui = -j: j (6) 

OUi 



and the dissipative phase fields dynamics follows from 

These equations lead in the limit e — > to the correct 
sharp interface limit above. For the case of static elastic- 
ity, this was carefully proved in 0. 

For the numerical realization, we employ explicit rep- 
resentations of both the elastodynamic equations and the 
phase field dynamics. The elastic displacements are de- 
fined on a staggered grid j3| ■ The derivation of the elas- 
todynamic equations of motion from a discretized action 
integral obeying invariance against time inversion guar- 
antees long time stability. 

We study crack growth in a rectangular geometry of 
a strip with fixed displacements at its upper and lower 
boundary. Horizontally, the grid can be shifted in order 
to keep the crack tip always in the center of the system. 
Thus, crack growth can be studied over long times in 
relatively small systems. Typical dimensions used here 
are 600 x 200 grid points, the phase field interface width 
is e = 5 Aa; (Ax is the lattice unit) and the Poisson ratio 
is v = 1/3. The Rayleigh speed vr is normalized to one. 
In the soft phase, we typically set the elastic constants 
to one percent of the values in the hard phase; however, 
these values are qualitatively not significant. Notice that 
after rescaling the equations of motion depend only on 
the driving force A and the kinetic coefficient D; in the 
numerical realization, also the phase field width e and 
the strip width L appear. 

First, we studied the growth of cracks in the vicinity of 
the Griffith point. Here, the tip radius in not determined 
by the length scale D/vr but by the phase field interface 
width e. In the strip geometry, the dimensionless driving 
force is A = u§(A + 2fi)/4Lj with the fixed vertical dis- 
placement u Q applied to the strip. As a nontrivial test, 
the numerical results validate the analytical prediction of 
the Griffith point A = 1 in the framework of the model, 
as the propagation velocity tends to zero (see Fig. [5J . 

The main goal was to approve that elastodynamics al- 
lows steady state growth without collapsing into the fi- 
nite time cusp singularity of the ATG instability, select- 
ing both a non-zero tip radius and a propagation velocity 
below the Rayleigh speed. The simulations confirm this 
prediction, and a typical steady state shape is shown in 
Fig-H Obviously, the tip radius is not determined by the 
intrinsic phase field length scale e. This can also be seen 
in Fig. |31 where we plotted the steady state tip radius r 
as function of the kinetic coefficient D for various driv- 
ing forces A. Only for very low kinetic coefficients, the 
tip radius is cut off by the interface width e, otherwise 
it is fairly bigger; for high kinetic coefficients the satura- 
tion is induced by the system size. In between, however, 
the scales are well separated and the radius ro is a lin- 
ear function of D, in good agreement with the theoreti- 
cal analysis: since both parameters v/vr and vro/D are 
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FIG. 2: Steady state velocity versus dimensionless driving 
force A; A = 1 is the Griffith point. The tip radius is 
here determined by the phase field interface width. We used 
D/tVR = 1.85 here. For A ^ 3.4 we get into the tip splitting 
regime. 
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FIG. 3: Tip radius ro as function of the kinetic coefficient 
D for different driving forces A. In the intermediate linear 
regime the length scales are well separated. 



FIG. 4: The quantity vro/D as function of the kinetic coeffi- 
cient for different driving forces A. 
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FIG. 5: vro/D as function of the driving forces A; the curves 
for different kinetic coefficients do not differ significantly, in 
agreement with the theoretical prediction. 
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t = 24.41 



predicted to be universal functions of the driving force 
A alone, we conclude that ro should depend linearly on 
the kinetic coefficient. Notice that also the tail opening 
A is of the scale A ~ D/vr instead of A ~ e. Further- 
more, the "dissipation rate" vro/D is rather independent 
of the kinetic coefficient (Fig. H). Thus the universal de- 
pendence of vro / D as a function of the driving force A 
can be extracted (Fig. |SJ. We believe that the results 
can be improved by increasing the system size and by a 
better separation of the length scales, which will be done 
in the near future. 

Snapshots of a typical tip splitting scenario for rela- 
tively high driving forces are shown in Fig. Repeated 
irregular splitting of the crack tip occurs, followed by 
symmetrical growth of the sidebranches. After a while, 
one finger wins the competition, moves back to the center 
of the strip and can finally split again. 

In summary, a phase field model has been developed 
to describe crack growth, based on thermodynamically 
well defined equations with a valid sharp interface limit. 
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FIG. 6: Irregular tip splitting scenario. We used D/eva = 
1.85 and A = 3.6; the system size is 600 x 400. Time is given 
in units D /v\. 
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The model shows the possibility of steady state growth 
of cracks and a tip splitting instability, in agreement with 
analytical predictions. In contrast to other models pre- 
viously discussed it provides a selection of the tip radius 
by the scale D/vr. 
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